Measuring the neutron star equation of state with gravitational wave observations 
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We report the results of a first study that uses numerical simulations to estimate the accuracy with which one 
can use gravitational wave observations of double neutron star inspiral to measure parameters of the neutron-star 
equation of state. The simulations use the evolution and initial-data codes of Shibata and Uryii to compute the 
last several orbits and the merger of neutron stars, with matter described by a parametrized equation of state. 
Previous work suggested the use of an effective cutoff frequency to place constraints on the equation of state. 
We find, however, that greater accuracy is obtained by measuring departures from the point-particle limit of the 
gravitational waveform produced during the late inspiral. As the stars approach their final plunge and merger, the 
gravitational wave phase accumulates more rapidly for smaller values of the neutron star compactness (the ratio 
of the mass of the neutron star to its radius). We estimate that realistic equations of state will lead to gravitational 
waveforms that are distinguishable from point particle inspirals at an effective distance (the distance to an 
optimally oriented and located system that would produce an equivalent waveform amplitude) of 100 Mpc or 
less. As Lattimer and Prakash observed, neutron-star radius is closely tied to the pressure at density not far 
above nuclear. Our results suggest that broadband gravitational wave observations at frequencies between 500 
and 1000 Hz will constrain this pressure, and we estimate the accuracy with which it can be measured. Related 
first estimates of radius measurability show that the radius can be determined to an accuracy of SR ~ 1 km at 
100 Mpc. 

PACS numbers: 04.25.dk, 04.30.Tv, 04.25.Nx, 26.60.Kp, 04.80.Nn 



I. INTRODUCTION 



Gravitational wave observations can potentially measure 
properties of neutron star equations of state (EOS) by mea- 
suring departures from the point-particle approximation to 
the gravitational waveform produced during the late inspiral 
phase of binary neutron star coalescence. We examine here 
the accuracy with which detectors with the sensitivity of Ad- 
vanced LIGO can extract from inspiral waveforms an EOS pa- 
rameter associated with the stiffness of the neutron star EOS 
above nuclear density. 

To do this, we use a set of numerical simulation waveforms 
produced by varying the EOS used to model the neutron star 
matter. The signal analysis focuses on the late inspiral, as the 
radius of the orbit r approaches the neutron star radius R. The 
orbital dynamics in this region will depend on the radius and 
internal structure of the neutron star, which in turn depend on 
the EOS. We also estimate the accuracy with which neutron 
star radii, closely linked to the EOS parameter varied, can be 
extracted. This is a preliminary study, using a first set of multi- 
orbit binary neutron star waveforms. Subsequent work will 
use an improved AMR code with higher resolution to obtain 
higher accuracy and to more fully explore the EOS parameter 
space. 

The study of radius and EOS effects on gravitational wave 
inspiral was first aimed at questions of detectability 0] [2j, 
showing that the tidal effects would only affect phase evolu- 
tion at the end of the inspiral and that point particle waveforms 



could be used for template-based detection in LIGO 1 . Subse- 
quently, the competing effects of relativity and finite size at 
the end of binary neutron star inspiral have been studied with 
an array of different approximations, yielding estimates of the 
gravitational wave spectrum which depart from that of a point 
particle inspiral starting somewhere between 500 Hz and over 
1000 Hz 015] ©Eli El Q31 

The study of EOS signature on gravitational waveforms has 
often focused on the merger and coalescence phases of the 
waveform lfl2l [131 H4l . above the frequency of the last or- 
bit. These frequencies are higher than the sensitive band of 
ground-based interferometers like Advanced LIGO, except in 
carefully tuned high frequency narrow-band configurations. 
Early work on measurability of finite size effects [15] used a 
model of point-particle inspiral truncated at such a frequency, 
which could be sought with narrow-band tuning. 

The characteristic frequency where EOS effects become 
important for the gravitational waveform is often estimated 
via an innermost stable circular orbit or Roche lobe overflow 
in a quasiequilibrium approximation. With many such esti- 
mates above 1000 Hz ||7] QT| [16), some summaries of neu- 
tron star radius measurability with gravitational waves have 
assumed that EOS dependence in double neutron-star binaries 
is unlikely to be detected with Advanced LIGO ifTTl . 

However, for quantitatively studying the late inspiral and 
merger phases of binary neutron stars, numerical relativity is 
required. Until quite recently, there has been no general rel- 



However with the increased sensitivity of Advanced LIGO, the contribu- 
tion to phase evolution from tidal effects may affect the waveform during 
the early, low-frequency part of the inspiral 1 3 1. 
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ativistic simulation that quantitatively clarifies the inspiral to 
merger phases primarily because of limitation of the computa- 
tional resources, although a number of simulations have been 
done for a qualitative study ll^[l3l[T8l[T^l2^l2Tll22ll23ll24l. 
The crucial drawback in the previous works was that the sim- 
ulations were short-term for the inspiral phase; the inspiral 
motion of the neutron stars was followed only for about one 
orbit, and the gravitational wave spectrum determined only 
above 1 kHz (but see Il25ll26l ). Furthermore, the simulations 
were usually started with a quasi-circular state in which the 
approaching velocity between two neutron stars is assumed to 
be zero and thus the eccentricity is not exactly zero. Thus, in 
the previous studies, the radial velocity at the onset of merger 
is not correctly taken into account and the non-zero eccentric- 
ity could play an unfavorable role. 

In this paper, by examining longer evolutions of binary neu- 
tron stars, the early relaxation from the quasi-circular state can 
be removed. By comparing the numerical inspiral waveform 
to that of point particles, we confirm that the frequency evo- 
lution in the late inspiral differs from the point particle case, 
accumulating phase more quickly in the final orbits and merg- 
ing at earlier times. This effect depends systematically on the 
EOS and resultant radius of the neutron stars simulated, and 
for large variations in EOS and radius the effect is larger than 
estimates of error in the numerical waveform. 

We calculate the signal strength of this difference in wave- 
form using the sensitivity curves of commissioned and pro- 
posed gravitational wave detectors, and find that there is a 
measurably different signal at reasonable distances from the 
inspirals of binary neutron stars with different EOS. This 
leads to a first quantitative estimate of the measurability of 
EOS with Advanced LIGO. 

We also note that the broadband configuration of Advanced 
LIGO does as well or better as narrow band configurations in 
detecting the difference between neutron star EOS. An im- 
proved understanding of the sensitivity of different gravita- 
tional wave detector configurations to neutron star structure 
will be essential for the design of next-generation detectors 
for gravitational wave astrophysics. 



n. EQUATIONS OF STATE 

We choose EOS based on work done in E71 to develop a 
parameterized EOS that accurately reproduces features of re- 
alistic EOS. Systematic variation of the EOS parameters al- 
lows us to determine which properties significantly affect the 
gravitational radiation produced, and thus can be constrained 
with sufficiently strong gravitational wave detections. 

The models chosen for this study use a variation of one 
EOS parameter in the neutron star core. The EOS pressure 
p is specified as a function of rest mass density 2 . Rest mass 



FIG. 1 : Initial choices of EOS for numerical evolution compared to 
the set of tabled EOS considered in 1271 . Candidates are labelled in 
order of increasing softness: 2H, H, HB, B, 2B. 
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density p, and the resulting energy density e are then deter- 
mined by the first law of thermodynamics. We use piecewise 
poly tropic EOS, of the form 



p(p) = K lP Y 



d- = —pd~, 
P P 



(1) 



Rest mass density p = m^n is proportional to number density n with the 
mass per baryon, if the matter were to be dispersed to infinity, fixed to be 
m b = 1.66 x 10~ 24 g. 



in a set of intervals pi < p < pi+x in rest mass density, with 
e/p — ► 1 as p — ► 0. 

A fixed crust EOS models the behavior of matter down to 
10 12 gem -3 ; the numerical simulations considered do not re- 
solve densities below this. The crust is modelled with a single 
polytrope region, fitted to tabulated crust EOS, for the region 
above neutron drip. The polytrope has r crust = 1.3569, with 
K crast chosen sothatp/c 2 = 1.5689 x 10 31 dyncm~ 2 when 
p = 10 13 gem -3 . The core EOS is constructed independently 
of crust behavior, and the dividing density between the crust 
and core varies by EOS. 

In 11271 it was found that three zones within the core are 
needed to accurately model the full range of proposed EOS 
models; however in this paper we will consider only one core 
zone, described by just one poly tropic EOS. We vary the core 
EOS with an overall pressure shift px, specified at the fidu- 
cial density p\ — 5.0119 x 10 14 gcm -3 , while holding the 
adiabatic index in all regions of the core fixed at T = 3. 
While only a subset of realistic EOS are well-approximated 
by a single core polytrope, reducing the EOS considered to 
this single-parameter family allows us to estimate parameter 
measurability with a reasonable number of simulations. 

After fixing the core adiabatic index, increasing the overall 
pressure scale pi produces a family of neutron stars with pro- 
gressively increasing radius for a given mass; the p\ parameter 
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TABLE I: Properties of initial EOS. These range from the "softest" 
EOS at the top, which results in a prompt collapse to a black hole 
upon merger, to the "hardest" (or "stiffest") at the bottom. Model 
HB is considered a typical EOS. The pressure pi, which is the pres- 
sure at density p\ = 5 x 10 14 g cm -3 , determines the polytropic EOS 
for the neutron star core; all candidates have Y = 3. Radius R and 
compactness GM/c 2 R are those of a single isolated M = 1.35 Mq 
TOV star where radius is measured in Schwarzschild-like coordi- 
nates. An astrophysically important EOS-dependent parameter is the 
maximum neutron star mass, M max , which is given in the fifth col- 
umn. 



Model log 10 pi [dyncm" 2 ] R [km] GM/c 2 R M max [M ] 



2H 


34.90 


15.2 


0.13 


2.83 


H 


34.50 


12.3 


0.16 


2.25 


HB 


34.40 


11.6 


0.17 


2.12 


B 


34.30 


10.9 


0.18 


2.00 


2B 


34.10 


9.7 


0.21 


1.78 



A. Initial data 

Conformally fiat initial data is generated by constructing a 
quasiequilibrium sequence of irrotational neutron stars in bi- 
nary system, following the methods of iFTl fT6l [3T1 [32l [33l [341 . 
As in previous work, we assume irrotational flow fields, ne- 
glecting spin of the neutron stars. This assumption is based 
on the estimation of negligible tidal spin-ups |fl~||2]]. The pa- 
rameterized EOS of Eq. ([T]l is incorporated in the code to solve 
for initial data with a conformally flat spatial geometry, using 
the Isenberg-Wilson-Mathews formulation 11351 [36l coupled 
to the neutron star matter equation consistently. Each of the 
stars has a baryon number equal to that of an isolated star with 
gravitational mass M — 1.35 Mq. The initial data for the full 
numerical simulation is taken from the quasi-equilibrium con- 
figuration at a separation such that ~ 3 orbits remain before 
merger. Relevant quantities of the initial configurations for 
each parameterized EOS are presented in Table [H] 



was chosen in part because Lattimer and Prakash [28 1 found 
that pressure at one to two times nuclear density is closely tied 
to neutron star radius, with R oc p^ 4 . The radius is less sensi- 
tive to variation of the adiabatic index in the neutron star core, 
for reasonable adiabatic indices |27l . 

An important element of future work will be incorporating 
additional variations of the EOS within the core. This could 
involve additional models of varying adiabatic index around a 
fixed pi, as well as multiple piecewise-polytrope zones within 
the core or EOS parameters yielding neutron stars of the same 
R but different internal structure. Such work would yield in- 
sight into the relative size and correlation of effects on the or- 
bital evolution due to the stellar radius and internal structure. 

The first models were chosen with EOS that "bracket" the 
range of existing candidates, seen in Fig. [T] These mod- 
els are HB with pi — 10 3440 dyncm~ 2 , a standard EOS; 
2H with Pl = 10 34 - 90 dyncm- 2 , a stiff EOS; 2B, with 
Pl /c 2 = 10 34 10 dyncm- 2 , a soft EOS. Additional mod- 
els B with pi/c 2 = 10 34 - 30 dyncm~ 2 and H with pi/c 2 — 
10 34 50 dyncm -2 , were chosen with small shifts in parameter 
from HB to better estimate local parameter dependence of the 
waveform. 



III. NUMERICAL METHODS 

For each EOS considered, we simulate the late inspiral and 
merger of a binary neutron star system. For this study, we 
fix the gravitational mass of each neutron star in the binary 
to 1.35M Q , an average value for pulsars observed in binary 
systems [29, 30]. We expect the significance of tidal effects 
in this configuration to be fairly representative of tidal effects 
over the relatively narrow range of masses and mass ratios 
expected in astrophysical binary neutron star system. 



B. Numerical evolution 

The Einstein equations are evolved with the original ver- 
sion of the Baumgarte-Shapiro-Shibata-Nakamura formula- 
tion ll37l [38ll in which we evolve the conformal factor, ip — 
(ln7)/12, the trace K of the extrinsic curvature, the confor- 
mal three metrics, 7^ = 7~ 1 / 3 7r ) -, the tracefree part of the ex- 
trinsic curvature, Aij = 7~ 1 / 3 (ify — Kjij/3), and an auxil- 
iary three-vector, Fj = S^djjik- Here 7^ is the three metric, 
Kij the extrinsic curvature, 7 = det(7y ), and K = Kirf 3 . 
As in 1 39], we evolve the conformal factor ip, not the inverse 
of ijj, because the cell-centered grid is adopted in our code, 
and hence, the black hole spacetime is handled in the mov- 
ing puncture framework |40l |4T1 . For the conditions on the 
lapse, a, and the shift vector, /3\ we adopt a dynamical gauge 
condition as in l39l . 

The numerical scheme for solving the Einstein equation 
is the same as that in ll39l : We use the fourth-order finite 
difference scheme in the spatial direction and a third-order 
Runge-Kutta scheme in the time integration, where the ad- 
vection terms such as (3 l di(p are evaluated by a fourth-order 
non-centered difference. 

The hydrodynamics equations are solved as in [139 j: We 
evolve p* = polu 1 ^ , Hi = hui, and e* = pau 1 — P/(pau t ), 
where p is the rest-mass density, ui is the three-component of 
the four velocity, P is the pressure, and h is the specific en- 
thalpy defined by h = 1+e+P/ p and e is the specific internal 
energy defined by e = e/p — 1. To handle advection terms in 
the hydrodynamic equations, a high-resolution central scheme 
B2ll is adopted with a third-order piecewise parabolic interpo- 
lation and with a steep min-mod limiter. In the present work, 
the limiter parameter, b, is set to be 2.5 (see (43 ] for detail 
about the parameter 6). 

The fluid evolution during inspiral is essentially free of 
shocks, and when there are no shocks the simulations use the 
cold EOS specified in Sec. [Il] During merger, when the evo- 
lution has shocks, we include a hot component with a thermal 
effective adiabatic index T\, as described in l43l . Shock heat- 
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TABLE II: Quantities of each initial data set for irrotational binary neutron stars. Each star has a baryon number Mo equal to that of an isolated 
M — 1.35 Mq star. The ADM mass Madm of the initial slice includes the binding energy, J is the total angular momentum of the initial slice. 
The binary compactness Co is defined by Co = (QMadmGc~ 3 ) 2 ^ 3 . 



Model 


Pmax [g Cm 3 ] 


M [Mq] 


M ADM [Mq] 


cJ/{GMl UM 


) fl/2-K [Hz] 


C 




2H 


3.73196 x 10 14 


1.45488 


2.67262 


0.993319 


324.704 


8.96966 x 10" 


-2 


H 


7.02661 x 10 14 


1.48385 


2.67080 


0.989524 


321.468 


8.90593 x 10" 


-2 


HB 


8.27673 x 10 14 


1.49273 


2.67290 


0.995361 


309.928 


8.69582 x 10" 


-2 


B 


9.77811 x 10 14 


1.50247 


2.67290 


0.992638 


314.170 


8.77522 x 10" 


-2 


2B 


1.38300 x 10 15 


1.52509 


2.67229 


0.987681 


321.170 


8.90375 x 10" 


-2 



ing in merger can increase the thermal energy up to ~ 20-30% 
of the total energy lfT2ll . 

Gravitational radiation is extracted both by spatially de- 
composing the metric perturbation about flat spacetime in the 
wave-zone with spin-2 weighted spherical harmonics and by 
calculating the outgoing part of the Weyl scalar ^4. For 
equal mass neutron stars (such as those studied here), the 
quadrupole (£ = 2,m = ±2) mode is much larger than any 
other mode; we consider just this mode in this analysis. 

The waveforms output from the simulations are the cross 
and plus amplitudes h + c 2 D/GM lot and h x c 2 D/GM t0( of the 
quadrupole waveform, as would be measured at large dis- 
tance D 3> GM tot /c 2 along the z axis perpendicular to the 
plane of the orbit, versus the retarded time t ret . Here M tot 
is the sum of the two neutron star masses when they are far 
apart, M tot = 2.7M Q . The strain is sampled at discrete val- 
ues evenly spaced in t, with a sampling rate At of between 
0.006 ms (for 2B) and 0.031 ms (for 2H) which depends on 
the time of simulation. There is junk radiation in the early part 
of the extracted waveforms for t mt ^ 0. We discard this part 
in the data analysis. 

Each simulation is typically performed for three grid reso- 
lutions. In the best-resolution case, the diameter of neutron 
stars is covered by 60 grid points. Convergence tests with dif- 
ferent grid resolutions indicate that, with the best grid resolu- 
tion, the time duration in the inspiral phase is underestimated 
by about 1 orbit. This is primarily due to the fact that angular 
momentum is spuriously lost by numerical dissipation. Thus, 
the inspiral gravitational waves include a phase error, and as a 
result, the amplitude of the spectrum for the inspiral phase is 
slightly underestimated. However, we find that the waveforms 
and resulting power spectrum for the late inspiral and merger 
phases, which we are most interested in for the present work, 
depend weakly on the grid resolution. 



IV. WAVEFORM ANALYSIS 



quency / of the quadrupole waveform is then estimated by 



We construct the complex quantity 

h = h-L. — ih-x 



(2) 



from the quadrupole waveform data. The amplitude and phase 
of this quantity define the instantaneous amplitude \h\ and 
phase <f> = ar g h of the waveform. The instantaneous fre- 



1 Acj> 
2tt A* 



(3) 



The numerical data can be shifted in phase and time by adding 
a time shift r to the time series points and multiplying the 
complex h by e 1 ^ to shift the overall wave phase by <j>. 

It is useful to define a reference time marking the end of the 
inspiral and onset of merger. A natural choice for the end of 
the inspiral portion, considering the behaviour of a PP inspiral 
waveform, is the time of the peak in the waveform amplitude 
\h\. However, the amplitude of the numerical waveforms os- 
cillates over the course of an orbit. A moving average of the 
waveform amplitude over 0.5 ms segments is used to average 
this oscillation; the end of the inspiral is then defined as the 
time at the end of the maximum amplitude interval. The re- 
sulting merger time i M will be marked by solid vertical lines 
in the plots to follow. 

The numerical waveforms begin at different orbital fre- 
quencies. To align them for comparison, they are each 
matched in the early inspiral region to the same post- 
Newtonian point-particle (or PP) waveform. 



A. Post-Newtonian point particle 

In full GR, point-particle inspiral is not well-defined, and 
one is left with the post-Newtonian point-particle (PP) ap- 
proximation and fully general relativistic black-hole numer- 
ical solutions as natural substitutes. Fortunately, the Taylor 
T4 3.0/3.5 post-Newtonian specification, introduced in PHI , 
agrees closely with numerical binary black hole waveforms 
for many cycles, up to and including the cycle before merger 
(see also B31 ). This empirical agreement allows us to adopt 
the Taylor T4 waveform as an appropriate PP baseline wave- 
form, compatible with full GR until the last cycles. We will 
show that the binary neutron star waveforms depart from this 
waveform 4-8 cycles (200-560 M tot ) before the best-fit PP 
merger. 

The TaylorT4 waveform is constructed by numerically in- 
tegrating for a post-Newtonian parameter x, which is related 
to the orbital frequency observed at infinity Q.: 



( GM tot n 

{ c 3 



2/3 



(4) 
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FIG. 2: Solid lines show numerical waveforms, scaled by 
c 2 D/GM tot , and aligned in time and phase to the same point-particle 
post-Newtonian inspiral (dashed line), using the method described in 
Sec. |IVB| The two dashed vertical bars indicate the portion of the 
waveform used for matching; the last vertical bar indicates the end 
of inspiral time tu for the numerical waveform. The top four simu- 
lations, 2H, H, HB, and B, show the start of post-merger oscillations 
from a hypermassive NS remnant in the simulation. 2B shows quasi- 
normal ringdown from a prompt collapse to a black hole following 
merger. 
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Following [44], we use 



dx 16 c 3 x 5 
dt ~YGMm 



487 



274 229 



' ' ^ 2 ' 72 576 



254 

■ir" 1 

1712 

~T05 



,5/2 



178 384 023 737 1475 



856 



33 530 011200 192 
3 310 



(5) 



•7- 



105 v ; ) 189 



7/2 



to find x(t) and then orbital phase of the binary is found from 

~3 



d<S> 
~dt 



GAL 



r 3/2 



(6) 



This yields a result to 3.5 post-Newtonian order for the phase 
evolution. The constants of integration are fixed by specifying 
the coalescence time i PP , when x — ► oo, and the orbital phase 
at this time <i>(i PP ) = $ PP , 

The amplitude of the quadrupole waveform is calculated 
to 3.0 order in [46 1. For the complex waveform hpp mea- 
sured along the orbital axis of the binary, using the appropriate 
spherical harmonic conventions, the result is 



G Mi 



D 



1 



-X 



j^r + 6^ 2 
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(7) 



where D is the distance to the observer, and logarithmic terms 
with a frequency scale have been absorbed into the the phase. 
The coalescence phase of the quadrupole point particle wave- 
form, </> PP = arg /i PP (i PP ), is determined by a choice of the 
orbital coalescence phase $ PP . 



B. Match to post-Newtonian point particle 

To match the numerical data to the PP inspiral waveform, a 
relative time shift and a relative phase shift must be specified 
by varying these parameters to obtain the best match. The 
masses of the point particles in the PP waveform are fixed to 
be the same as the neutron stars in the numerical simulations 
- the gravitational mass of isolated (TOV) neutron stars with 
the same number of baryons - and so masses are not varied 
in finding the best match. With a goal of signal analysis, we 
choose the time and phase shift by maximizing a correlation- 
based match between two waveforms. 

The complex numerical relativity derived quadrupole wave- 
form, /inr = /i+ R — i/ix R , is convolved with the complex post- 
Newtonian quadrupole waveform, hpp = h ? J—ih ? ^ of Eq. (|7j, 
over a given matching region T\ < t < Tp. 



where M tot is the sum of the point masses. To first post- 
Newtonian order, x ~ M tol /r where r is the orbital radius. 



z(t;Ti,T f )= / h NR (t)h FP (t - r)dt. 



(8) 
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FIG. 3: Phase-optimized, time-limited match between numerical 
inspiral waveforms and point particle post-Newtonian waveforms. 
Contours are shown at 0.95, 0.98, 0.99, and 0.997. as a function 
of match region truncation at some time before numerical merger, 
Tf—tu, and relative shift between numerical and point particle wave- 
forms, t^ F — tu. The start of the match is fixed to 1.5 ms after the 
start of the numerical waveform. Subsequent analysis in this paper 
is done using the best match at a fixed Tp — tu of 1.8 ms for each 
waveform. 
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This quantity is similar to the complex matched filter output 
of the FINDCHIRP algorithm 07). The real part of z(t; T h Tp) 
corresponds to the correlation of the two waveforms as a func- 
tion of the time shift r. The absolute value \z(t; T\, Tp) | is the 
correlation maximized over a constant overall phase shift, and 
the argument is the overall phase shift required to so maxi- 
mize. A normalized match between two waveforms is 



m(r;T I ,T F ) = 



z(r\T h T ¥ ) 



c N r(0; Ti, T F ) ct pp (t; T h T F ) 
where the normalizing constants are defined by 



a 2 (r;T h T P ) = f *\h{t-r)\ 



dt 



(9) 



(10) 



where h is either /i NR , which defines er NR (r; T\, Tp), or /i PP , 
which defines <t pp (t; Ti, Tp). 

The shift r maximizing the match between two waveforms 
is sensitive to the portion of the numerical waveform matched, 
Ti < t < Tp. Some truncation of the numerical waveform 
is required to eliminate residual effects of initial data. One 
would also like to truncate the waveform at some point before 
the peak amplitude, at the end of the region without signifi- 
cant finite size effects. To determine where this region is, we 
plot the match between numerical waveforms as a function of 
both the end point of the match Tp, and the shift r, measuring 
quantities relative to the previously defined merger time tu of 
each numerical waveform and the PP coalescence time t* p so 
thatr = t PP - t u , in Fig. |U 

Reassuringly, for most of the waveforms, once we truncate 
a millisecond or so before the numerical merger tu, there is 
a region where the range of well-matched shifts r show little 
dependence on the exact end point Tp. This continues with 



earlier end points until the segment of the waveform consid- 
ered becomes so short that the range of matching r broadens 
significantly. The exception is the waveform of HB, which 
was the first simulation of the series. HB exhibits somewhat 
larger eccentricity, and drift in the best match timeshift, com- 
pared to the other waveforms. 

For subsequent analysis, we choose to take r maximizing 
the match for a region of the numerical waveforms between 
1.5 ms after the start of the waveform and 1.8 ms before nu- 
merical merger. Varying the details of this choice can change 
the best match t by up to ~ 1 ms. By comparison, one wave- 
form cycle takes between 0.5 and 2 ms in the inspiral region, 
so this will be a significant source of uncertainty in SNR es- 
timates. Longer or more accurate simulations are required to 
more precisely fix a post-Newtonian match. 



C. Comparison of waveforms 

Unlike the case of matching binary black hole simulations 
to point particle post-Newtonian ll44l |45|, the binary neutron 
star simulations show departure from point particle many cy- 
cles before the post-Newtonian merger time. Fig.|2]shows the 
four waveforms shifted so the best-match PP waveforms have 
the same i PP and cf> PP '. As the stiffness of the EOS and thus 
the radius of the neutron stars, increases, the end of inspiral 
for the binary neutron stars is shifted away from the end of 
inspiral for point particle post-Newtonian. 

This can also be seen by plotting the instantaneous fre- 
quency of the numerical simulation waveform with the same 
time shifts, as in Fig. [4] which also shows more clearly the 
difference in the post-merger oscillation frequencies of the 
hyper-massive remnants, when present. The larger neutron 
star produced by the stiff EOS 2H has a lower oscillation fre- 
quency than that from the medium EOS HB. The remnant 
forms with a bar-mode oscillation stable over a longer pe- 
riod than the ~ 10 ms simulated. The signal from such a bar 
mode may be even stronger when the full lifetime is included, 
although aspects of the physics neglected in this study will 
likely come into play. Information that can be extracted from 
the presence (or absence) and characteristics of a post-merger 
oscillation signal would complement the information present 
in the late inspiral. This is another subject for future study. 



D. Frequency spectrum of waveforms 

Given h+ or h x , one can construct the discrete Fourier 
transforms (DFTs) h + or h x . Both polarizations yield the 
same DFT amplitude spectrum \h\, with phase shifted by 7r/2, 
if one neglects discretization, windowing, and numerical ef- 
fects (including eccentricity). The amplitude spectrum \h\ is 
independent of phase and time shifts of the waveform. 

The stationary phase approximation is valid for the post- 
Newtonian waveform up to frequencies of about 1500 Hz 
(with < 10% error), so is used to plot the amplitude of the 
point particle spectrum. In terms of an instantaneous fre- 
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FIG. 4: Time-frequency behavior, vertical line markings as previous 
figure. The departure from the point particle time-frequency rela- 
tions, shown using a long-dashed line, occurs between 700-1000 Hz 
depending on EOS. 
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(11) 



the Fourier transform of the waveform has the amplitude 



\h\ ~ A(f) 



-1/2 



(12) 



For binaries comprised of equal mass companions, the gravi- 
tational radiation is dominated by quadrupole modes through- 
out the inspiral. The wave phase is negligibly different from 
twice the orbital phase 2<I> until the onset of merger at very 
high frequencies so we can use the relation 



df 
dt 



1 



7T GM 2 



,1/2 



dx 
~dt'' 



(13) 



to write the amplitude of the Fourier transform entirely in 
terms of the functions dx/dt of Eq. pi and the amplitude of 
the polarization waveforms A(f) — \h\ in Eq. (j7|. 

The translation of emitted waveforms into the strain ampli- 
tude measured at a detector involves transformations incorpo- 
rating the effects of the emitting binary's angle of inclination 
and sky location. These effects are absorbed into an effective 
distance D e g, which is equal to the actual distance D for a bi- 
nary with optimal orientation and sky location, and is greater 
than the actual distance for a system that is not optimally ori- 
ented or located. The detector will detect a single polariza- 
tion of the waveform, some combination of the plus and cross 
polarizations of the emitted waveform. The polarizations ex- 
tracted from the simulation can be used as two estimates of the 
strain at the detector for a given numerically modelled source, 
which give very close results and are subsequently averaged. 

To compare with noise curves in the usual units, the quan- 
tity / 1//,2 |/i(/)| is plotted, at a reference distance of D e ff = 
100 Mpc and rescaling from previously plotted numerical out- 
put h(t)c 2 D/GM tot using M tot = 2.7M Q . 

The full spectra of models 2H, H, HB, and B, seen in Fig. [5] 
show peaks at post-merger oscillation frequencies; those of H 
and B are weaker as the waveform is truncated shortly after the 
formation of the hyper-massive remnant. Waveforms of 2H 
and HB are also truncated while the post-merger oscillation 
is ongoing; if the simulations were allowed to continue, these 
peaks would presumably grow further. The simulation of 2B, 
in contrast, collapses to a black hole and has only a short lived 
(and relatively small-amplitude) quasinormal mode ringdown 
post merger. 

Note that time-frequency plots like the ones shown Fig. [4] 
showed that numerical waveforms follow the PP waveform at 
instantaneous frequencies of up to 700-1000 Hz, depending 
on the EOS. The disagreement in the spectra in Fig. [5] from 
the PP stationary-phase approximation waveform at frequen- 
cies below this is primarily due to the finite starting time of the 
numerical waveforms. To estimate spectra from the full inspi- 
ral, we construct hybrid waveforms. The short-term numer- 
ical waveforms are are smoothly merged on to long-inspiral 
PP-PN waveforms with Hann windowing to smoothly turn on 



w(n) 



1 — cos 



N —1 



(14) 
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FIG. 5: DFT of full numerical waveforms, at an effective distance D c ff = 100 Mpc, compared to noise spectra for Advanced LIGO (labelled 
"AdvLIGO" for the standard configuration and "Broadband" for the broad-band configuration) and the Einstein Telescope (labelled "ET") 
shown by thick grey lines. The DFT of the numerical waveforms turned off after £m is shown by dot-dashed lines, the stationary-phase point 
particle is shown by a dashed line for reference. The lower right figure shows a combined plot of inspiral-truncated waveforms, smoothly 
joined on to best-match PP inspiral time series before the DFT is taken. 
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or turn off 



w(n) = 



1 



1 



cos 



N — 1 



(15) 



a signal over a range of N points, < n < N (or over a time 
N At). To construct hybrid waveforms, these windows are 
used to turn on the numerical waveform as the matched post- 
Newtonian is turned off, such that the sum of the two window 
functions is 1 over the matching range. 

Although the post-merger oscillations are an interesting 
source of potentially measurable strains, the dependence on 
the cold EOS is less straightforward as temperature effects be- 
come significant during the the merger. We focus instead on 
the signal from the waveforms during the inspiral region only, 
turning off the waveforms after the end of inspiral, tyi- The 
DFT amplitudes of inspiral-only hybrid waveforms are also 
plotted in Fig. [5] from the range of EOS considered. We see 
here that HB is estimated to depart from PP earlier than either 



H or B, rather than at the expected intermediate value. This is 
ascribed to the higher eccentricity and lower accuracy of the 
match for the HB waveform. 



V. DETECTABILITY 

The question for neutron star astrophysics is whether these 
differences in the waveform will be measurable. We consider 
the possibility of detecting EOS effects with Advanced LIGO 
style detectors in varying configurations. 

We use several detector configurations commonly con- 
sidered for Advanced LIGO include tunings optimized for 
1.4 Mq NS-NS inspiral detection ("Standard"), for burst de- 
tection ("Broadband"), and for pulsars at 1150 Hz ("Narrow- 
band"). The sensitivity is expressed in terms of the one-sided 
strain-equivalent amplitude spectral density Sh(f) (which has 
units of Hz -1 / 2 ) of the instrumental noise in Advanced LIGO. 
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We also consider a provisional noise curve for the Einstein 
Telescope [4-8 1 . We consider only a single detector of each 
type, rather than a combination of detectors, for a preliminary 
estimate of detectability. 

Given two signals hi and hi, define the usual inner product 
for a given noise spectrum Sh(f) ll49l : 



(hi\ha) =4Re 



hi(f)K(f) 
Sh(f) 



df. 



(16) 



This inner product yields a natural metric on a space of wave- 
forms with distance between waveforms weighted by the in- 
verse of the noise. We filter the detector output s against an 
expected waveform h using (s\h). Then 



(17) 



is the optimal statistic to detect a waveform of known form 
h in the signal s. If the detector output contains a particular 
waveform h that is exactly matched by the template used, then 
the expectation value of g is the expected signal-to-noise ratio 
or SNR of that waveform: 



Q : 



(18) 



Given a signal hi to be used as a template, one can consider 
whether a known departure from this signal can be measured. 
Assuming the modified waveform h% is known, the expected 
SNR of /12 — hi is similarly 



£?diff 



= \Z(k) - hi\h 2 - hi). 



(19) 



The two signals will be (marginally) distinguishable 3 , if the 
difference between the waveforms has TJdiff > 1. 

We analyse the measurability of the differences between 
hybrid inspiral-only waveforms matched to point particle PN 
waveforms in the early inspiral. After the numerical wave- 
forms have been matched to the same post-Newtonian point- 
particle inspiral, the signals will be aligned in time and phase. 
We can then then compare the resulting waveforms to each 
other, and to a PN-only waveform, using different Advanced 
LIGO noise spectra. 

We report results in terms of the SNR measured by a single 
detector at an effective distance D e s = 100 Mpc. Results can 
be scaled to any distance; as each measured h oc X/D, Q&g oc 
1/D. 

The difference between waveforms due to finite size effects 
is not detectable in the NS-NS detection optimized configu- 
ration of Advanced LIGO for ~ 100 Mpc effective distances. 
However, in both narrowband and broadband the differences 
can be significant, and waveforms are distinguishable from 
each other and from the PP waveform. Note that this im- 
plies that even if the details of a NS-NS late inspiral signal 
are not known, the difference between it and a point particle 



TABLE III: £>diff in standard (NS-NS detection optimized) noise 
x (100Mpc/Aff) 

Model 2B B HB H 2H 



PP 0.32 0.45 0.55 0.46 0.69 

2B 0.36 0.48 0.38 0.63 

B 0.21 0.12 0.58 

HB 0.27 0.60 

H 0.58 



TABLE IV: 

x (100Mpc/Aff) 



in broadband (burst-optimized) noise 



Model 2B B HB H 2H 



PP 

2B 
B 

HB 
H 



1.86 2.32 2.67 2.38 2.89 
1.92 2.32 2.03 2.54 
0.81 0.80 2.27 







1.28 2.37 
2.35 



' Compare discussion in 1 50 1 of indistinguishability. 



waveform should be measurable. The quantity gdiff between 
the observed waveform and a best fit point particle waveform, 
limited to differences at high frequency, may be useful in it- 
self to constrain possible EOS without reference to waveform 
details. 



A. Parameter estimation 

We will assume that EOS effects on the waveform impact 
the late inspiral only. For simplicity, we assume that orbital 
parameters, such as M tot , mass ratio rj, point particle post- 
Newtonian coalescence time t£ p , and phase shift cA PP , are de- 
termined from the observations of the earlier inspiral wave- 
form, with sufficient accuracy that their measurement uncer- 
tainty will not affect the accuracy to which the late inspiral 
effects determine the EOS parameters. These measurements 
would be made by a broad-band instrument, in which the 
signal-to-noise ratio is expected to be high (~ 40 at 100 Mpc) 
and measurement accuracy is expected to be good 11491 . In- 
accuracies in these measurement could lead to biases in the 
measured EOS. This will be an important aspect to assess 
when high-quality binary neutron star simulations with vari- 
ous masses become abundant. 

With a one-parameter family of waveforms sampled, we 
can estimate the accuracy to which this parameter can be mea- 
sured. There are also other EOS-related parameters which are 
not considered. In the direct analysis of the measurability of 
the EOS parameter pi, we ignore variations of T within the 
core. Correspondingly, in the analysis of radius measurement, 
variations of the internal structures are neglected. Expanding 
coverage of the EOS parameter space is an ongoing project. 
However, these initial parameter choices are expected to give 
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TABLE V: gdS in narrowband 1150Hz noise x (lOOMpc/A: 



Model 2B B HB H 2H 

PP 0.91 2.75 3.69 2.65 2.12 

2B 1.92 2.92 1.82 1.45 

B 1.14 0.22 1.43 

HB 1.34 2.25 

H 1.42 



TABLE VII: 51og(pi), where p\ is measured in dyncm~ 2 , in broad- 
band (burst-optimized) noise x (D c ff/100Mpc). Compare the range 
in pi of the candidates from Tableland Fig.[T] 



Model 2B 


B H 2H 


2B - 


0.10 0.20 0.32 


B 


- 0.25 0.26 


H 


- 0.17 



TABLE VI: SR (km) in broadband (burst-optimized) noise 
x (LW/lOOMpc) 



Model 2B B H 2H 



2B 
B 
H 



0.63 1.28 2.17 
- 1.74 1.89 
- 1.23 



the dominant contributions to finite size effects of the wave- 
form 4 . 

We estimate errors in parameter estimation to first order in 
1/g or, equivalently, in S6 A , using the Fisher matrix Tab — 
(d A h\d B h) ED- Its inverse, (T- 1 )^, yields 



-1\AB 



S0 A S9 B = (T- 1 ) 
so that the expected error in a given parameter 6 A is 



(S9 A y = (r- 1 ) 



1\AA 



(20) 



(21) 



and the cross terms of the inverse Fisher matrix yield correla- 
tions between different parameters. 

With a few simulations of varying parameter value, we es- 
timate dh/dpi and dh/dR from two of the sampled wave- 
forms, hi and h 2 . For a single parameter 9 (which can be 
taken to be either pi or R), we have 



dh 
09 



~ — 

8=(8 1 +0 2 )/2 ^2 



hi 



(22) 



where hi = h(9i) and hi = h(6 2 ), and then, for our one- 
parameter family where we neglect correlations with other pa- 
rameters, we have to first order 



(S9Y 



(92 - 9if 



(h 2 - hi\h 2 - hi) ' 



(23) 



Using adjacent pairs of models to estimate waveform de- 
pendence at an average parameter value, we then find esti- 



mates of radius measurability as shown in Table VI and pi 



4 For example, using the 1PN tidal effect estimates of | 3 51] yields about 
10% variation in the tidal terms contributing to binding energy and lumi- 



nosity from changing internal structure- 
while keeping radius fixed. 



-varying the apsidal constant- 



measurability as shown in Table |VII| for the burst-optimized 
noise configuration. 



B. Some sources of error 

We have cavalierly neglected many higher order but likely 
relevant effects in this preliminary analysis. 

It is possible that tidal effects measurably influence the or- 
bital evolution before the start of the numerical simulations, 
as estimated in 0, slowly enough not to be seen over the few 
cycles of the waveform matched to PP in this analysis. In one 
sense this analysis is a worst-case scenario, as it assumes exact 
PP behavior before the numerical match. Earlier drift away 
from point particle dynamics would give larger differences 
between waveforms, and more sensitive radius measurement, 
but poses a challenge by requiring accurate numerical simula- 
tion over many cycles to verify EOS effects. Combining nu- 
merical estimation with PN analyses like those of [3] and/or 
quasiequilibrium sequence information may clarify the transi- 
tion between effectively PP and tidally influenced regimes. 

The waveforms themselves have some residual eccentricity 
from initial data and finite numerical resolution. One can esti- 
mate eccentricity error by comparing the plus polarization to 
the cross polarization shifted by 7r/2 from the same numerical 
waveform. This results in f^ffi of ~ 0.3 for HB and 2H, rather 
than the expected quadrupole polarization cross-correlation of 
zero. (The other waveforms produce ^Jiff ^0.05 between po- 
larizations). 

The length of the inspirals, as discussed in the section on 
PN matching, limits precision in choosing the best match 
time. We can estimate these effects on current results by vary- 
ing the match region considered; this changes 'Q&ff results by 
up to ~ 0.5 at 100 Mpc in the broadband detector. The reso- 
lution from existing numerical simulations is thus comparable 
to the difference between parameters of the three closest mod- 
els. 

We have only a coarsely sampled family of waveforms; es- 
timates of dh/d9 are limited by this. The value of q^s for 
HB-B and H-HB should be half that of H-B, but instead they 
are the same or greater — we are hitting the limit of numeri- 
cal and matching accuracy. We can also estimate the validity 
of the linear parameter dependence in the central models by 
comparing HB to (H + B) /2. This results in JSff — 0-8, an- 
other estimate of systematic error. 

We conclude that, although these are of course first esti- 
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mates, they should be better than order-of-magnitude; the un- 
certainty in each g^g is comparable to the q&b between the 
H/HB/B waveforms, and smaller than the ~g^ff for larger dif- 
ferences in EOS and for comparison with PP. 

Finally, we note that use of a Fisher matrix estimate of pa- 
rameter measurement accuracy is fully valid only in high SNR 
limit of ftJiff > 10 ||52ll , i.e., for distances < 20 Mpc in the 
broadband detector. The results do not take into account mul- 
tiple detectors, nor multiple observations, nor parameter cor- 
relations. A full estimation of EOS parameter measurability 
will require more detailed analysis, with a larger set of inspi- 
ral simulations sampling a broader region of parameter space, 
e.g., with mass ratios departing from unity, and with several 
more orbits prior to merger. 

VI. CONCLUSIONS 

Gravitational wave astrophysics provides a promising new 
window on behavior of cold dense matter. We estimate that 
realistic EOS will lead to gravitational inspiral waveforms 
which are distinguishable from point particle inspirals at an 
effective distance of 100 Mpc or less in a burst-optimized Ad- 
vanced LIGO configuration, as good or better than in narrow 
band detector configuration. 

While the standard noise configuration of Advanced LIGO 
is not sensitive to the differences in the waveform, the prelim- 
inary standard noise curve of the Einstein Telescope indicate 
the ability to differentiate between EOS at roughly double the 



distance as broadband Advanced LIGO. In general, detun- 
ing detectors to be more sensitive at frequencies above 700 Hz 
will lead to improved gravitational wave constraints on neu- 
tron star EOS and radius. 

First estimates of parameter measurability in broadband 
Advanced LIGO show 5R ~ 1km x (100 Mpc/L» eff ). One 
can also consider a direct constraint on the EOS pressure pa- 
rameter pi at a rest mass density p\ = 5 x 10 14 gcm -3 
of Spi ~ 10 32 dyncm -2 at an effective distance D e s — 
100 Mpc. These estimates neglect correlations between these 
parameters and other details of internal structure, but such de- 
tails are expected to give relatively small corrections to the 
tidal effects. 

While our results must still be considered preliminary, they 
strongly motivate further work on gravitational wave con- 
straints from binary neutron star inspirals. Future numerical 
simulations with longer inspirals and increased coverage of 
parameter space should improve the accuracy of the estimates. 
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